Navier-Stokes with a body and varying free stream



In [1]:

    
using ViscousFlow



In [2]:

    
using Plots
pyplot()
default(grid = false)

Solve flow past a stationary body

Set the flow parameters



In [3]:

    
Re = 1000; # Reynolds number

Set the free stream motion



In [4]:

    
U = 1.0; # Mean free stream velocity
σx = 0.4; # amplitude of velocity fluctuation 
ϕ = -π/2  # phase lag of pitch to heave
fstar = 0.25/π # fc/U

K = π*fstar # reduced frequency, K = πfc/U

oscil = RigidBodyTools.OscilX(2π*fstar,U,U*σx,ϕ);
U∞ = RigidBodyMotion(oscil)









    Out[4]:





Rigid Body Motion:
  ċ = 0.6 + 0.0im
  c̈ = 0.0 + 0.0im
  α̇ = 0.0
  α̈ = 0.0
  RigidBodyTools.OscilX(0.5, 1.0, 0.4, -1.5707963267948966, Constant (0.0), AddedProfiles:
  Constant (1.0)
  0.4 × (Sinusoid (ω = 0.5) >> 3.141592653589793)
, d/dt (AddedProfiles:
  Constant (1.0)
  0.4 × (Sinusoid (ω = 0.5) >> 3.141592653589793)
))

Set up points on the body. Here is a plate:



In [5]:

    
n = 100;
body = Plate(1.0,n)









    Out[5]:





Plate with 100 points and length 1.0 and thickness 0.0
   Current position: (0.0,0.0)
   Current angle (rad): 0.0



In [6]:

    
cent = 1.0 + 1.0im
α₀ = -20*π/180 # angle of attack
T = RigidTransform(cent,α₀)
T(body) # transform the body to the current configuration









    Out[6]:





Plate with 100 points and length 1.0 and thickness 0.0
   Current position: (1.0,1.0)
   Current angle (rad): -0.3490658503988659

Set up the domain



In [7]:

    
xlim = (0.0,4.0)
ylim = (0.0,2.0)









    Out[7]:





(0.0, 2.0)



In [8]:

    
plot(body,xlim=xlim,ylim=ylim)









    Out[8]:

Now set up the coordinate data for operator construction



In [9]:

    
X = VectorData(body.x,body.y);
X̃ = VectorData(body.x̃,body.ỹ);

Set the domain size and time step size



In [11]:

    
Δx = 0.01;
Δt = min(0.5*Δx,0.5*Δx^2*Re);

Now set up the system

Set up the state vector and constraint force vector for a static body



In [14]:

    
sys = NavierStokes(Re,Δx,xlim,ylim,Δt, X̃ = X, isstore = true)









    Out[14]:





Navier-Stokes system on a grid of size 416 x 208



In [15]:

    
w₀ = Nodes(Dual,size(sys));
f = VectorData(X);
xg, yg = coordinates(w₀,dx=Δx,I0=origin(sys))









    Out[15]:





(-0.075:0.01:4.075, -0.035:0.01:2.035)

Set up the integrator here



In [16]:

    
plan_intfact(t,u) = CartesianGrids.plan_intfact(t,u,sys)
plan_constraints(u,t) = ConstrainedSystems.plan_constraints(u,t,sys)
r₁(u,t) = ConstrainedSystems.r₁(u,t,sys,U∞)
r₂(u,t) = ConstrainedSystems.r₂(u,t,sys,U∞)

@time solver = IFHERK(w₀,f,sys.Δt,plan_intfact,plan_constraints,(r₁,r₂),rk=ConstrainedSystems.RK31)









    



 15.137296 seconds (158.01 M allocations: 8.069 GiB, 9.45% gc time)






    Out[16]:





Order-3 IF-HERK integrator with
   State of type Nodes{Dual,416,208,Float64,Array{Float64,2}}
   Force of type VectorData{100,Float64,Array{Float64,1}}
   Time step size 0.005

Initialize the state vector and the history vectors



In [17]:

    
t = 0.0
w₀ .= 0.0
u = deepcopy(w₀)

fx = Float64[];
fy = Float64[];
thist = Float64[];

uhist = [];
tsample = 0.2; # rate at which to store field data

Advance the system!

Set the time range to integrate over.



In [18]:

    
tf = 1.0;
T = Δt:Δt:tf;



In [19]:

    
for ti in T
    global t, u, f = solver(t,u)
    
    push!(thist,t)
    push!(fx,sum(f.u)*Δx^2)
    push!(fy,sum(f.v)*Δx^2)
    (isapprox(mod(t,tsample),0,atol=1e-6) || isapprox(mod(t,tsample),tsample,atol=1e-6)) ? push!(uhist,deepcopy(u)) : nothing
end
println("solution completed through time t = ",t)









    



solution completed through time t = 1.0000000000000007

Plotting

Basic plot



In [20]:

    
plot(xg,yg,vorticity(uhist[end],sys),levels=range(-15,15,length=30), color = :RdBu,clim=(-15,15))
plot!(body)









    Out[20]:

Make a movie



In [24]:

    
@gif for i = 1:length(uhist)
    plot(xg,yg,vorticity(uhist[i],sys),levels=range(-15,15,length=30), color = :RdBu,clim=(-15,15))
    plot!(body)
end









    












    



┌ Info: Saved animation to 
│   fn = /Users/jeff/JuliaProjects/ViscousFlow/examples/tmp.gif
└ @ Plots /Users/jeff/.julia/packages/Plots/sbXPh/src/animation.jl:104






    Out[24]:



In [23]:

    
px = plot(thist,2*fx,xlim=(0,Inf),ylim=(0,3),xlabel="Convective time",ylabel="\$C_D\$",legend=false)
py = plot(thist,2*fy,xlim=(0,Inf),ylim=(-3,3),xlabel="Convective time",ylabel="\$C_L\$",legend=false)
plot(px,py)









    Out[23]: